\chapter{Introduction}
\label{intro} 
In this thesis we present methods for biophysically constrained, image based estimation of myocardial motion.
 Myocardial motion is a primary measure of myocardial function, and differences in myocardial function can help diagnose
 and quantify cardiomyopathies better than myocardial structural differences.  The main feature of our method
 is that we estimate the active contraction of the myocardial fibers instead of the displacements. This is accomplished with
 a biophysically-constrained, four-dimensional (space plus time) formulation that aims to complement information that can
 be gathered from the images by a priori knowledge of cardiac mechanics and electrophysiology.

\section{Motivation}

Cardiac diseases claim more lives than any other disease in the world \cite{aha}. Early diagnosis and
treatment can save many lives and reduce the associated socio-economic costs. Cardiac diseases
are characterized by both changes in the myocardial structure as well as changes in cardiac
function. Consequently, it is important for any cardiac diagnosis method to consider both these
aspects. 

Imaging can help in the diagnosis of cardiac masses, cardiomyopathy, myocardial infarcion, and valvular disease.  Advances in Cine Magnetic Resonance (MR) imaging methods have enabled us to acquire high-resolution 4D images of the heart that capture the structural and functional characteristics of individual hearts. Studies suggest that standardizing acquisition protocols and objective analysis especially of regional myocardial function will help improve the accuracy and reduce inter-observer variability \cite{bluemke03}. This has
created the need for sophisticated and highly automated image analysis methods, which can
identify and precisely quantify subtle and spatially complex patterns of functional
changes in the heart. The development of such methods is the primary goal of this thesis.


% Alternative specialized pulse
% sequences (e.g., tagged cine-MRI) can create rich image features, but
% they have lower resolution and signal-to-noise ratio. In addition,
% computational challenges limit analysis to 3D (2D $\times$ time)
% motion estimation, where in fact 4D analysis would be preferable
% \cite{odonell-e06,sampath-prince07}. 
% 
% An alternate approach is to estimate myocardial motion from MR Cine sequences. MR Cine images in a clinical setting at sub millimeter resolutions (in-plane), with slice thickness in the range of 6-10mm. The temporal resolution varies between 25-70ms. A lot of work has been done in extracting cardiac motion fields from MR Cine sequences [15, 22-27]. 

% These can be classified into two main categories. The first approach uses segmentation of the myocardial wall, followed by geometric and mechanical modeling using active contours or surfaces to extract the displacement field and to perform the motion analysis [15, 24, 27]. For matching two contours or surfaces, curvatures are frequently used to establish initial sparse correspondences, followed by the dense correspondence interpolation in other myocardial positions by regularization or mechanical modeling [15, 23]. The lack of distinct landmarks on the myocardial wall makes it difficult to estimate the wall motion based on surface tracking. In addition this approach is very sensitive to the accuracy with which the myocardium can be segmented. Also it performs poorly in regions within the myocardium, and manages to only align the myocardial boundaries. The other approach uses energy-based warping or optical flow techniques to compute the displacement of the myocardium [22, 25, 26]. Perperidis et al. [25] use a regular grid with a B-spline basis to model the deformation and use normalized mutual information as the similarity metric which is calculated over the whole image. One of the major shortcomings of these approaches is that the transformation estimated as a result of the registration is not unique and in fact does not necessarily conform to the underlying myocardial motion. The same algorithm can give different estimates of motion for different initial guesses and different parameters. The problem arises since these methods attempt to maximize the image similarity with only a smoothness constraint on the transformation. Since there can be many transformation that can map an image onto another (especially sparsely sampled ones as in the case of MR Cines) there is no guarantee that the estimated transformation is the correct one [28, 29]. These methods estimate motion by evolving the current estimate of motion under the action of external image forces (image similarity) and internal forces which constrain the regularity of the motion (smoothness). Such regularizers work well with respect to noise removal but they do not incorporate a priori knowledge of the underlying cardiac motion. Incorporating a biomechanically-inspired model for the myocardium has the potential for a more accurate motion estimation [30]. Functional models of the heart are direct computational models, designed to reproduce in a realistic manner the cardiac activity, often requiring high computational costs and the manual tuning of a very large set of parameters. Such methods can be computationally prohibitive for our purposes, and we instead select a level of modeling compatible with reasonable computing times and involving a limited number of parameters. Such simplifications add additional modeling errors, but our hypothesis is that in spite of these modeling errors the estimated motion fields shall be more accurate than those obtained from approaches not incorporating a priori knowledge. A detailed and thorough review of cardiac image registration methods can be found in [31] and a general review of image registration methods can be found in [32]. 

% Segmentation of the ventricles and the myocardium is the first step
% toward quantitative functional analysis of cine-MRI data.  However,
% segmentation is time consuming, thereby limiting clinical throughput
% \cite{axel02}. Moreover, sometimes accuracy is limited by long-axis
% motion, and inter and intra-observer variability\cite{bluemke03}.

\section{Overview}

Cardiomyopathies are characterized both by structural changes in the myocardial structure as well as changes in myocardial function. It is important to consider both these changes when designing an automatic diagnosis algorithm. The method by which the structural and functional information is obtained from the patient should be simple and account for minimal patient discomfort and trauma. Magnetic resonance imaging (MRI) has been very popular in the recent years for acquiring structural information in the human body. Apart from providing us with structural information of the myocardium, these images can provide us with myocardial wall motion if they can be extracted accurately. 

The task of extracting myocardial wall motion from MR Cine sequences is not trivial, and although a lot of work has been done in this field, most of the approaches rely purely on image similarity to drive the motion estimation and therefore fail to successfully capture the true underlying cardiac motion. As is common in image registration, regularization (typically $L_2$) is usually used to constrain the motion estimation process. When used for motion estimation, only a smoothing regularization like $L_2$ is not sufficient and we are not guaranteed to get the true underlying motion. This is illustrated in Figure \ref{fig:twist}, where we see that two different motions can produce the same configuration. Because of this reason, it is important to provide intelligent priors to help the motion estimation algorithm obtain physically correct motions. In this thesis we present a framework for incorporating physics based constraints, in the form of partial differential equations (pde), for the motion estimation problem. 



% To handle the computational issues with such models, we present algorithms for parallel algorithms for solving pde systems on adaptive octree meshes.     
  
\begin{figure}
	
	\subfloat[] {  
		\includegraphics{new_tikz/twist1}
		}
	\subfloat[] {  
		\includegraphics{new_tikz/twist2}
	}	
	\subfloat[] {  
		\includegraphics{new_tikz/twist3}
	}
\caption{A schematic of the left ventricle is shown in (a). The same contracted configuration can be reached via two deformations; (b) via radial contractions, and (c) via a twist. }	
\label{fig:twist}		
\end{figure}



\section{Related work}

Cardiomyopathies present themselves in different forms, both by structural changes like plaque formation, fat deposits, etc., and functional changes like variations in ventricular wall motion, ejection fraction, and perfusion. Both of these need to be extracted from the image before accurate diagnosis can be done. A lot of work has been done in feature extractors for structural characterizations of disease. This has focused primarily on extracting image features like intensities and gradients \cite{intgrad}, moments \cite{hammer,moments}, Gabor features \cite{manju96}, and local frequency representations \cite{locfreq}. The problem with cardiomyopathies is that not all of them can be characterized by structural changes. Function at rest may be abnormal as a result of one of the spectrum of ischemic heart diseases (ischemia, infarction, hibernation) or of cardiomyopathy from other causes. During stress testing, new or worsening wall motion abnormalities are indicative of functionally significant coronary artery stenosis \cite{smart00}. In addition, wall motion imaging to detect regional contractile reserve is an accurate measure of myocardial viability, and the results can help guide coronary revascularization therapy. Characterizing cardiomyopathies based on both structural and functional changes will make the diagnosis algorithm more accurate and robust. Of course quantization of the myocardial wall motion represents a challenge in itself. 

Most clinical modalities used to image myocardial function evaluate passive wall motion (ventriculography) or wall thickening (echocardiography, gated single-photon emission computed tomography, or cine MR imaging). MR imaging also allows quantitative measurement of regional intramyocardial motion and, subsequently, strain, which can be more sensitive to wall motion abnormalities than is wall thickening. MR imaging methods for the quantification of intramyocardial wall motion can be loosely classified into two approaches, those relying on specially developed MR imaging protocols to help in the estimation of myocardial motion and those relying on image analysis techniques to extract motion estimates from MR Cine sequences.

\subsection{Specialized MR Protocols}
\subsubsection{MR Tagging}
MR Tagging was developed to provide non-invasive virtual markers inside the myocardium, which deform with myocardial motion \cite{mrtag}. MR imaging and especially tagged MR are currently the reference modalities to estimate dense cardiac displacement fields with high spatial resolution. The deformation fields, as well as the derived motion parameters such as myocardial strain can be determined within an accuracy of 2mm x 2mm \cite{Shi99,chenBook}. The primary disadvantage of tagging is the reduced spatial resolution of strain relative to the image spatial resolution. In tagging, after the displaced tag lines are detected \cite{guttman94}, the displacement field can be estimated and intramyocardial strain can be computed in a variety of ways \cite{mrtag}. With this approach, although strain may be interpolated to any desired spatial resolution, the fundamental spatial resolution of strain is nominally determined by the distance between the tag lines, which is typically several pixels. Tag detection has an additional disadvantage in that it typically requires substantial manual intervention and is therefore a time-consuming task. Harmonic phase analysis will likely obviate tag detection \cite{osman99}, but the spatial resolution of the resultant strain maps will not necessarily improve. The spatial resolution of strain maps obtained from tagged images after harmonic phase analysis is determined by the k-space filter of the analysis; in practice with single breath-hold acquisitions, the resolution has been relatively poor \cite{garot00}. Additionally since the right ventricle (RV) is much thinner than the left ventricle (LV), it is difficult to place more than a single tag within the RV, making the estimation of RV motion extremely difficult and inaccurate. Since we are most interested in the characterizing RV function, tagging is not appropriate for our purpose.
\subsubsection{Phase Contrast Imaging}
The second approach is that of MR phase contrast imaging \cite{mrphase}, which is based on the concept that spins that are moving in the same direction as a magnetic field gradient develop a phase shift that is proportional to the velocity of the spins. This information can be used directly to determine the velocity of the spins, or in the cardiac case the velocity of any point within the myocardium. The main problem with this approach is that four acquisitions have to be made for each heart, one the regular MR cine sequence and one phase contrast acquisition each for the velocity components in the x, y, and the z directions. Consequently, MR phase contrast imaging is not used much in a clinical setting. 
\subsubsection{DENSE and HARP}
Displacement encoded imaging with stimulated echoes (DENSE) \cite{epstein04} and harmonic phase imaging (HARP) \cite{osman99} employ 1-1 spatial modulation of magnetization to cosine modulate the longitudinal magnetization as a function of position at end diastole. Later in the cardiac cycle the cosine-modulated signal is sampled and used to compute myocardial strain from the signal phase. The sampled signal generally includes three distinct echoes:  a displacement encoded stimulated echo, the complex conjugate of the displacement-encoded echo, and an echo arising from T1 relaxation. If the T1 relaxation and complex conjugate echoes are suppressed, then a phase image representing just the displacement encoded echo can be reconstructed. However, data acquisition in single breathhold DENSE MR imaging has been limited to only one cardiac phase. Multiple breathhold DENSE produces images at multiple phases of the cardiac cycle, but the resolution has been poor and is fundamentally a 2D approach and the estimation of through-plane displacement has been poor. 

\subsection{Extracting motion from MR Cine images}
An alternate approach is to estimate myocardial motion from MR Cine sequences. MR Cine images in a clinical setting at sub millimeter resolutions (in-plane), with slice thickness in the range of $6-10$mm. The temporal resolution varies between $25-70$ms. A lot of work has been done in extracting cardiac motion fields from MR Cine sequences \cite{ledesma01,McE00,papademetris-01,perperidis04,Shi99,Song91,Wang01}. These can be classified into two main categories. The first approach uses segmentation of the myocardial wall, followed by geometric and mechanical modeling using active contours or surfaces to extract the displacement field and to perform the motion analysis \cite{papademetris-01,Shi99,Wang01}. For matching two contours or surfaces, curvatures are frequently used to establish initial sparse correspondences, followed by the dense correspondence interpolation in other myocardial positions by regularization or mechanical modeling \cite{McE00,Shi99}. The lack of distinct landmarks on the myocardial wall makes it difficult to estimate the wall motion based on surface tracking. In addition this approach is very sensitive to the accuracy with which the myocardium can be segmented. Also it performs poorly in regions within the myocardium, and manages to only align the myocardial boundaries. The other approach uses energy-based warping or optical flow techniques to compute the displacement of the myocardium \cite{ledesma01,perperidis04,Song91}. Perperidis et al. \cite{perperidis04} use a regular grid with a B-spline basis to model the deformation and use normalized mutual information as the similarity metric which is calculated over the whole image. One of the major shortcomings of these approaches is that the transformation estimated as a result of the registration is not unique and in fact does not necessarily conform to the underlying myocardial motion. The same algorithm can give different estimates of motion for different initial guesses and different parameters. The problem arises since these methods attempt to maximize the image similarity with only a smoothness constraint on the transformation. Since there can be many transformation that can map an image onto another (especially sparsely sampled ones as in the case of MR Cines) there is no guarantee that the estimated transformation is the correct one \cite{Cachier:MICCAI:01,Cachier-JMIV-2004}. These methods estimate motion by evolving the current estimate of motion under the action of external image forces (image similarity) and internal forces which constrain the regularity of the motion (smoothness). Such regularizers work well with respect to noise removal but they do not incorporate a priori knowledge of the underlying cardiac motion. Incorporating a biomechanically-inspired model for the myocardium has the potential for a more accurate motion estimation \cite{mcculloch98}. Functional models of the heart are direct computational models, designed to reproduce in a realistic manner the cardiac activity, often requiring high computational costs and the manual tuning of a very large set of parameters. Such methods can be computationally prohibitive for our purposes, and we instead select a level of modeling compatible with reasonable computing times and involving a limited number of parameters. Such simplifications add additional modeling errors, but our hypothesis is that in spite of these modeling errors the estimated motion fields shall be more accurate than those obtained from approaches not incorporating a priori knowledge. A detailed and thorough review of cardiac image registration methods can be found in \cite{makela02} and a general review of image registration methods can be found in \cite{Zitova03}. 

\subsection{Motion estimation using biomechanical models}
 To address these problems in motion reconstruction, one of the main thrusts in recent research has been 4D
motion estimation using biomechanical models.
There is significant work on the integration of imaging with
cardiovascular mechanics. In \cite{hu-metaxas-axel-03}, a piecewise
linear composite biomechanical model was used to determine active
forces and the strains in the heart based on tagged MRI information.
In \cite{papademetris-01} and \cite{papademetris02}, echocardiography
and MR images were combined with biomechanical models for cardiac
motion estimation.  Interactive segmentation was combined with a
Bayesian estimation framework that was regularized by an anisotropic,
passive, linearly elastic, myocardium model.  The authors recognized
the importance of neglecting active contraction of the left ventricle.
In \cite{sermesant-02a,sermesant-e06,sermesant-e06a}, the need for
realistic simulations and the need for inversion and data assimilation
was outlined. In \cite{moireau-chapelle-letallec08}, Kalman filters
were used to recover the initial state of the heart and spatial
abnormalities. That method however, is difficult to generalize to
nonlinear inversion with time-dependent inversion parameters.


\section{Contributions}

The major contributions of this thesis are,

% \begin{enumerate}[i.]
\begin{itemize}
  \item Development of an automatic method for estimation patient specific myocyte orientations using template warping of {\em in vitro} diffusion tensor images.
  	\item A parallel bottom-up algorithm for coarsening octrees,
	which is also used for partitioning the input in our other
	algorithms.
	\item A parallel bottom-up algorithm for constructing linear
	octrees. We avoid the synchronization issues that are usually
	associated with parallel top-down methods.
	\item An algorithm for enforcing 2:1 balance refinement in
	parallel. The algorithm constructs the minimum number of nodes to
	satisfy the 2:1 constraint. Its key feature is that it avoids
	parallel searches, which as we show in sections \ref{sec:parBal}
	and \ref{sec:commCost}, are the main hurdles in achieving good
	isogranular scalability.
  \item  An algorithm to build element-to-nodal connectivity information efficiently and to compress this information in a way that achieves a three-fold compression (a total of  four words per octant). These methods are explained in sections \ref {sec:oda} and \ref{sec:compress}. 
	\item A lookup table based conforming discretization scheme that requires
only a single traversal for the evaluation of a partial differential
operator. 
	\item A novel and efficient way of computing mutual information. The proposed method defines a non-uniform, adaptive sampling scheme for estimating the entropies of the images, which is less vulnerable to local maxima as compared to uniform and random sampling.
	\item A patient-specific image-based inversion formulation for the active forces; 
	\item A multigrid-accelerated, octree-based, adaptive finite-element forward solver that incorporates anatomically-accurate fiber information; and 
	\item An adjoint/Hessian-based inversion algorithm. 
\end{itemize}
%\end{enumerate}

Our work on estimating patient specififc myocyte orientations was published in \cite{sundar06a}. The work on octree construction and 2:1 balance refinement and meshing in \cite{sundar-biros-05,sundar-sampath-biros-e07}. Our work on image similarity measures using octrees has been published in \cite{sundar07}. 

\section{Limitations}

The limitations of our current implementation (but not the method) are
the assumptions of linear geometric and material response and the
potential bias due to template-based fibers that does not account for
anatomical variability, that is still requires some preprocessing of
the initial frame to assign material properties and fiber orientation,
that assumes zero residual stresses and initial conditions, and that
it does not include an electrophysiology model.

Our on-going work includes transition to an intensity-based
image-registration inversion (in which case we need to solve a
nonlinear inversion).  Among the many open problems are the
level of required model complexity for clinically relevant motion
reconstructions, the bias of the fibers, the sensitivity to the values
of the material properties, and the sensitivity to the image
similarity functional.

\section{Organization of this thesis}

In the rest of the thesis, we discuss algorithms for bio-physically constrained cardiac motion estimation. We start by describing the bio-mechanical model of the heart and issues related to the estimation of myocyte orientations in Chapter \ref{chap:model}. In Chapter \ref{chap:octree} we present algorithms for the construction, 2:1 balance refinement, meshing and a PETSc \cite{petsc-web-page} compatible distributed array for solving linear systems. We also present our algorithm for estimating mutual information using octrees in this chapter. Finally, in Chapter \ref{chap:inverse} we present the pde constrained optimization framework for estimation of myocardial motion.

%Each chapter has an extensive discussion on the related work and technically precise discussion on the contributions.